% code used to create Fig. 9 in "An alternative pathway for delivery of 
% power by outer hair cells to cochlear traveling waves", Samaras and Meaud
% Plots the nonlinear response vs frequency of model 4 with revised MET conductance function
% the response is plotted at multiple location
% Julien Meaud
% julien.meaud@me.gatech.edu
% ------------------------------------------------------------------------

close all
clear all

%% Colors
orange = [0.8500, 0.3250, 0.0980];
green = [140, 170, 0]./255;
yellow  = [0.9290, 0.6940, 0.1250];
purple = [0.4940, 0.1840, 0.5560];
grey = [0.6 0.6 0.6];
black = [0 0 0];
magenta = [1 0 1];
blue_CB = [26 133 255] ./ 255;
purple_CB = [212 17 89] ./ 255;
panelLabels = {'A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P'};

plotRL=0;
plotNormReCF=0;

% Patches for shading
xPatch = [0 55 55 0];
yPatch3 = [-10 -10 0 0];
yPatch4 = [0 0 2 2];
transparency_Patch = 0.5;

plotBMblack=0;
if plotBMblack 
    BMcolor = [0 0 0];
else
    BMcolor=blue_CB;
end
LS1 = '-';
LW=1.5;

leftX=0.12;
gapX=0.02;
height=0.25;
width=0.2;

Level=[30,50,70,80];

nLevel = length(Level);
transp = linspace(0.3,1,nLevel);

deltaOHC = 0;
fSize=10;

plotNormModel=1e3;
if plotNormModel==1e3
    FreqLabel='Freq. (kHz)';     
else
    FreqLabel='Freq. re BF';
end

CFV=[2500,8000,20000,35000];  % CF of the plotted location
r = 3;  % number of rows (1st row: BM and OHC/DC amp; 2nd row: BM and TM amp; 3rd row: phase re BM)
c = length(CFV); ; % number of columns

xTickCF=[1000 2500 3500;4000 8000 11000; 10000 20000 25000;15000 35000 45000]/1000;
xTickCF=xTickCF';

for iLevel=1:length(Level)

    Level_curr=Level(iLevel);
       
    if iLevel==1
        freqV=[1000:500:1000 1500:250:3000 3500:500:8000 9000:1000:50000];
    elseif iLevel==2
        freqV=[1000:500:1000 1500:250:3000 3500:500:7000,8000:1000:50000];
    elseif iLevel==3
        freqV=[1000:500:1500 2000:250:3000 3500:500:6000,7000:1000:50000];
    elseif iLevel==4
        freqV=[1000:500:1500,2000:250:3000 4000:500:7000,8000:1000:50000];
    end
    
    % load the data
    clear  uBM_M uDC_M uTMB_M
    for iFreq=1:length(freqV)
        freq=freqV(iFreq);

        load(sprintf('TimeDomain_ModifiedM4/PTTD_G36_v13_%gPT%g_SMALL.mat',Level_curr,freq))              
                
        xMesh=ModelGeneratedVarsStruct.xMesh;
        freq=StimulusOptionsStruct.pureTone.F;
    
        uBM=ModelOutputsPostFFT.uBM;
        uDC=ModelOutputsPostFFT.uDC;  
        uTMB=ModelOutputsPostFFT.uTMB;
    
        uBM_M(:,:,iFreq)=uBM;
        uDC_M(:,:,iFreq)=uDC;
        uTMB_M(:,:,iFreq)=uTMB;
    
    end

        % find the best place using the low level data
    
  for iCF=1:length(CFV)
    CF=CFV(iCF);
    
    if iLevel==1
        indexFreq=find(freqV==CF);    
        [~,iBP(iCF)]=max(abs(uBM_M(:,2,indexFreq)));  
    end

    uBMreECP=squeeze(uBM_M(iBP(iCF),2,:)*1e7./convertSPLtoPressureInPascals(Level_curr)).*cos(deltaOHC);
    uDCreECP=squeeze(uDC_M(iBP(iCF),2,:)*1e7./convertSPLtoPressureInPascals(Level_curr));
    uTMBreECP=squeeze(uTMB_M(iBP(iCF),2,:)*1e7./convertSPLtoPressureInPascals(Level_curr));
        
        fig1=figure(1);
        set(fig1,'units','inches','position',[2 2 7 5.5])
        subplot('Position',[leftX+(width+gapX)*(iCF-1) 0.71 width height])
        plot(freqV./plotNormModel,20*log10(abs(uDCreECP)),'Color',[magenta transp(iLevel)],'LineWidth',LW,'Linestyle',LS1,'handlevisibility','off')
        hold on
        plot(freqV./plotNormModel,20*log10(abs(uBMreECP)),'Color',[BMcolor transp(iLevel)],'LineWidth',LW,'Linestyle',LS1,'Displayname',sprintf('%g dB',Level_curr))
        
        subplot('Position',[leftX+(width+gapX)*(iCF-1) 0.4 width height])
        plot(freqV./plotNormModel,20*log10(abs(squeeze(uBMreECP))),'Color',[BMcolor transp(iLevel)],'LineWidth',LW,'Linestyle',LS1,'Displayname',sprintf('%g dB',Level_curr))
        hold on
        plot(freqV./plotNormModel,20*log10(abs(squeeze(uTMBreECP))),'Color',[green transp(iLevel)],'LineWidth',LW,'Linestyle',LS1)
        
        subplot('Position',[leftX+(width+gapX)*(iCF-1) 0.1 width height])
        plot(freqV./plotNormModel,unwrap(angle(squeeze(uDC_M(iBP(iCF),2,:)./uBM_M(iBP(iCF),2,:))))/(2*pi),'Color',[magenta transp(iLevel)],'LineWidth',LW)
        hold on
        plot(freqV./plotNormModel,unwrap(angle(squeeze(uTMB_M(iBP(iCF),2,:)./uBM_M(iBP(iCF),2,:))))/(2*pi),'Color',[green transp(iLevel)],'LineWidth',LW,'Linestyle',LS1)

  end
end
    
 
% process the passive data (linear simulations) ---------------------------
for iCF=1:length(CFV)
     CF=CFV(iCF);
    if iCF==1      
        vars_Pass = load('Models_M1M2M3M4_FreqDomain/PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_5_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat');     
    end

    uBM_Pas = vars_Pass.ModelOutputsStruct.structural.uBM(iBP(iCF),:);
    uTMB_Pas = vars_Pass.ModelOutputsStruct.structural.uTMB(iBP(iCF),:);
    uTMR_Pas = vars_Pass.ModelOutputsStruct.structural.uTMR(iBP(iCF),:);
    freqPas = vars_Pass.StimulusOptionsStruct.freq;
   
    uBM_PasReECP = uBM_Pas.*1e7.*cos(deltaOHC)./convertSPLtoPressureInPascals(20);   
    uTM_PasReECP = uTMB_Pas.*1e7./convertSPLtoPressureInPascals(20);  
    uDC_Pas = vars_Pass.ModelOutputsStruct.structural.uDC(iBP(iCF),:);
    uDC_PasReECP = uDC_Pas.*1e7./convertSPLtoPressureInPascals(20);
    
    figure(1)
    subplot('Position',[leftX+(width+gapX)*(iCF-1) 0.71 width height])
    plot(freqPas./plotNormModel,20*log10(abs(uBM_PasReECP)),'Color',[BMcolor],'LineWidth',LW,'Linestyle',':','Displayname','Passive')
    hold on
    if iCF==1
        h2 = legend('AutoUpdate','off','Orientation','horizontal');
        h2.Location = 'northwest';
     end
     plot(freqPas./plotNormModel,20*log10(abs(uDC_PasReECP)),'Color',[magenta transp(iLevel)],'LineWidth',LW,'Linestyle',':')               
    
    if iCF==1      
         text(1.08*CF/1000,57,'BM','Color',BMcolor,'Fontsize',10)
         text(0.25*CF/1000,20,'OHC-DC','Color',magenta,'Fontsize',10)
    end

    xlim([5 27.5]*1e3./plotNormModel*CF/20000)
    ylim([0 90])
    yticks([0:20:80])
    yticklabels([0:20:80])
    if ~plotNormReCF
        xticks( xTickCF(:,iCF))
        xticklabels([])
    else
        xticks( xTickCF(:,iCF))
        xticklabels([]);
    end
    text(24000./plotNormModel*CF/20000,83,panelLabels{iCF},'Fontsize',16,'Fontweight','bold')
    yticks([0:20:80])
    if iCF== 1
        ylabel(['Disp. re ECP',char(10),'(dB re 1 nm/Pa)'])
        yticklabels([0:20:80])
       
    else
        yticklabels([])
    end

    vline(CF./plotNormModel,'k--')
    set(gca,'Fontsize',fSize)
    title(['x=',num2str(xMesh(iBP(iCF))),' cm'])
       
    subplot('Position',[leftX+(width+gapX)*(iCF-1) 0.40 width height])
    plot(freqPas./plotNormModel,20*log10(abs(uBM_PasReECP)),'Color',[BMcolor],'LineWidth',LW,'Linestyle',':','Displayname','Passive')
    hold on
   
    if plotRL
        plot(freqPas./plotNormModel,20*log10(abs(uRL_PasReECP)),'Color',[purple_CB transp(iLevel)],'LineWidth',LW,'Linestyle',':')
    else % plot TM
        plot(freqPas./plotNormModel,20*log10(abs(uTM_PasReECP)),'Color',[green transp(iLevel)],'LineWidth',LW,'Linestyle',':')
    end
    xlim([5 27.5]*1e3./plotNormModel*CF/20000)
    ylim([0 90])
    yticks([0:20:80])
    yticklabels([0:20:80])
    if ~plotNormReCF
       xticks( xTickCF(:,iCF))
       xticklabels([])
    else
       xticks([0.25:0.25:1.75])
       xticklabels([]);
    end
    text(24000./plotNormModel*CF/20000,83,panelLabels{iCF+c},'Fontsize',15,'Fontweight','bold')
    yticks([0:20:80])
    if iCF == 1
        ylabel(['Disp. re ECP',char(10),'(dB re 1 nm/Pa)'])
        yticklabels([0:20:80])
    else
        yticklabels([])
    end
    vline(CF./plotNormModel,'k--')
    set(gca,'Fontsize',fSize)
    if iCF==1
        text(7000./plotNormModel*CF/20000,50,'TM','Color',green,'Fontsize',10)
        text(15000./plotNormModel*CF/20000,22,'BM','Color',BMcolor,'Fontsize',10)
    end
    
    subplot('Position',[leftX+(width+gapX)*(iCF-1) 0.1 width height])
    plot(freqPas./plotNormModel,unwrap(angle(uDC_PasReECP./uBM_PasReECP))./(2*pi),'Color',[magenta transp(iLevel)],'LineWidth',LW,'Linestyle',':')
    hold on
    plot(freqPas./plotNormModel,unwrap(angle(uTM_PasReECP./uBM_PasReECP))./(2*pi),'Color',[green  transp(iLevel)],'LineWidth',LW,'Linestyle',':')    
    xlim([5 27.5]*1e3./plotNormModel*CF/20000)
    ylim([-0.5 0.5])
      
    yticks([-1:0.25:1])
    if iCF==1
        yticklabels([-1:0.25:1])
    else           
      yticklabels([])
    end
        
    if ~plotNormReCF
        xticks( xTickCF(:,iCF))
        xticklabels( xTickCF(:,iCF))
    else
        xticks([0.25:0.25:1.75])
        xticklabels([0.25:0.25:1.75])            
    end
    
    text(24000./plotNormModel*CF/20000,0.4,panelLabels{iCF+2*c},'Fontsize',16,'Fontweight','bold')
    vline(CF./plotNormModel,'k--')
    set(gca,'Fontsize',fSize)

    xlabel(FreqLabel)
    if iCF==1
        ylabel(['Phase re BM',char(10),'(cycles)'])
    end
        
end 
